## calculate term structure of SDRs based on estimates of UC model for real rate

rm(list=ls())
graphics.off()
set.seed(616)
library(dplyr)
source("R/sdr_fns.R")

## CHOOSE HERE
## spec <- "y1"        # one-year real rate
spec <- "y10"    # ten-year real rate

## settings
LB <- 0
year1 <- 1990
year2 <- 2019

## load estimation results
filename <- paste0("results/uc_estimates_", spec, ".RData")
load(filename)
print(range(data$year))

t1 <- which(data$year == year1)
t2 <- which(data$year == year2)
tbl <- cbind(data$rr[c(t1,t2)], data$rstar.mean[c(t1,t2)])
tbl <- cbind(tbl, tbl[,1] - tbl[,2])
colnames(tbl) <- c("Real rate", "r*", "r-tilde")
rownames(tbl) <- c(year1, year2)
print(round(tbl, 2))
rstar_new <- mean(X[,1,t2])
rstar_old <- mean(X[,1,t1])

## term structure of discount rates
N <- 400 # maximum maturity
mats <- 1:N             # yield maturities in years (for plotting)
cat("Simulating from posterior distribution...\n")
P1 <- sim_uc_sdr(X[,1,t1], sqrt(sigeta2), sqrt(sigv2), phi, LB=LB)
P2 <- sim_uc_sdr(X[,1,t2], sqrt(sigeta2), sqrt(sigv2), phi, LB=LB)
y1 <- -100/(1:N)*log(P1)
f1 <- -100*diff(log(c(1,P1)))
y2 <- -100/(1:N)*log(P2)
f2 <- -100*diff(log(c(1,P2)))  # or: f <- Psim[1:(N-1)]/Psim[2:N]-1

## figure with shifting term structures
dev.new()
par(mar=c(3,3,.5,.5), mgp=c(2,.6,0))
cols <- c("red", "steelblue")
plot(mats, y1, type="l", lwd=2, ylim=range(0,y1, y2), xlab="Horizon (years)", ylab="Discount rate (percent)", col=cols[1], main="")
lines(mats, y2, lwd=2, col=cols[2])
abline(h=rstar_old, lty=2, lwd=2, col=cols[1])
abline(h=rstar_new, lty=2, lwd=2, col=cols[2])
legend("topright", legend=c(year1, year2), lwd=2, col=cols, bg="white")

## save term structures of SDRs
save(file=paste0("results/sdr_uc_", spec, ".RData"), rstar_old, rstar_new, P1, f1, y1, P2, f2, y2, N)
